/*
 * GeneralParticleSource.cc
 *
 * @date Jan 18, 2011
 * @author Tim Niggemann, III Phys. Inst. A, RWTH Aachen University
 * @copyright GNU General Public License v3.0
 */

#include "ParticleSource.hh"

#include <G4ParticleTable.hh>
#include <G4UImanager.hh>
#include <G4RunManager.hh>
#include <G4Run.hh>
#include <G4OpticalPhoton.hh>
#include <Randomize.hh>

#include <CLHEP/Units/SystemOfUnits.h>

#include <boost/any.hpp>
#include <boost/format.hpp>
#include <fstream>
#include <algorithm>

#include "ParticleSourceMessenger.hh"
#include "rnd/Direction.hh"
#include "rnd/Flat.hh"
#include "rnd/General.hh"
#include "rnd/ThreeVector.hh"
#include "rnd/FixedValue.hh"

ParticleSource::ParticleSource()
{
	particleGun = new G4ParticleGun(1);
	messenger = ParticleSourceMessenger::getInstance();
}

ParticleSource::~ParticleSource()
{
	//
}
// 计算光子的偏振
G4ThreeVector ParticleSource::getPolarizationVector(double angle, G4ThreeVector kphoton)
{
	// Calculate polarization.
	G4ThreeVector normal(1., 0., 0.);
	G4ThreeVector product = normal.cross(kphoton);
	G4double modul2 = product * product;
	G4ThreeVector e_perpend(0., 0., 1.);
	if (modul2 > 0.)
	{
		e_perpend = (1. / sqrt(modul2)) * product;
	}
	G4ThreeVector e_paralle = e_perpend.cross(kphoton);
	G4ThreeVector polar = cos(angle) * e_paralle + sin(angle) * e_perpend;
	return polar;
}
// 将位置或动量方向从法线为(0,0,-1)的平面转换到源的表面平面
G4ThreeVector ParticleSource::transformToSourcePlane(G4ThreeVector v)
{
	const G4ThreeVector dicingPlaneNormal(0, 0, 1);
	const G4ThreeVector sourcePlaneNormal = messenger->getSurfaceNormal();
	const G4ThreeVector rotationAxis = dicingPlaneNormal.cross(sourcePlaneNormal);
	G4double rotationAngle = 0.;
	if (rotationAxis.mag() != 0.)
	{
		rotationAngle = acos(
			dicingPlaneNormal * sourcePlaneNormal / (dicingPlaneNormal.mag() * sourcePlaneNormal.mag()));
	}
	else if (sourcePlaneNormal == -dicingPlaneNormal)
	{
		rotationAngle = 180. * CLHEP::deg;
	}
	if (rotationAngle == 180. * CLHEP::deg)
	{
		v *= -1;
	}
	else if (rotationAngle != 0.)
	{
		v.rotate(rotationAxis, rotationAngle);
	}
	return v;
}

void ParticleSource::GeneratePrimaries(G4Event *event)
{
	// const unsigned int nParticles = messenger->getNParticles();
	rnd::RandomDouble *particles = (rnd::RandomDouble *)new rnd::General("resources/psf_number_distribution.properties");
	const unsigned int nParticles = (int)particles->shoot();
	// Init random number generators.
	// 从messenger中获取输入的位置分布，如果没有输入，则使用默认的Flat分布
	messenger->setA(1 * CLHEP::mm);
	messenger->setB(1 * CLHEP::mm);
	rnd::ThreeVector *position = new rnd::ThreeVector(new rnd::Flat(-messenger->getA() / 2., messenger->getA() / 2.),
													  new rnd::Flat(-messenger->getB() / 2., messenger->getB() / 2.), new rnd::FixedValue(0.0, 0.0));
	rnd::RandomDouble *polar = new rnd::FixedValue(messenger->getPolar(), 360.0 * CLHEP::deg);
	// 从messenger中获取输入的能量分布，如果没有输入，则使用默认的Flat分布
	// rnd::RandomDouble* energy =
	// 		messenger->getEInput().empty() ?
	// 				(rnd::RandomDouble*) new rnd::Flat(messenger->getEMin(), messenger->getEMax()) :
	// 				(rnd::RandomDouble*) new rnd::General(messenger->getEInput());
	// 从文件中读取能量分布
	rnd::RandomDouble *energy = (rnd::RandomDouble *)new rnd::General("resources/psf_energy_distribution.properties");
	// Update messenger values.
	// 更新messenger中参数
	if (!messenger->getEInput().empty())
	{
		G4UImanager::GetUIpointer()->ApplyCommand(
			boost::str(boost::format("/ps/energy/eMin %g eV") % (((rnd::General *)energy)->getXMin() / CLHEP::eV)));
		G4UImanager::GetUIpointer()->ApplyCommand(
			boost::str(boost::format("/ps/energy/eMax %g eV") % (((rnd::General *)energy)->getXMax() / CLHEP::eV)));
	}

	// 从messenger中获取输入的时间分布，如果没有输入，则使用默认的Flat分布
	// rnd::RandomDouble* time =
	// 		messenger->getTInput().empty() ?
	// 				(rnd::RandomDouble*) new rnd::Flat(messenger->getTMin(), messenger->getTMax()) :
	// 				(rnd::RandomDouble*) new rnd::General(messenger->getTInput());
	// 从文件中读取时间分布
	rnd::RandomDouble *time = (rnd::RandomDouble *)new rnd::General("resources/psf_time_distribution.properties");
	// Update messenger values.
	// 更新messenger中参数
	G4double psf_TMin = ((rnd::General *)time)->getXMin() / CLHEP::ns;
	G4double psf_TMax = ((rnd::General *)time)->getXMax() / CLHEP::ns;
	if (!messenger->getTInput().empty())
	{
		G4UImanager::GetUIpointer()->ApplyCommand(
			boost::str(boost::format("/ps/tMin %g ns") % (((rnd::General *)time)->getXMin() / CLHEP::ns)));
		G4UImanager::GetUIpointer()->ApplyCommand(
			boost::str(boost::format("/ps/tMax %g ns") % (((rnd::General *)time)->getXMax() / CLHEP::ns)));
	}
	G4UImanager::GetUIpointer()->ApplyCommand("/g4sipm/persist/tMin " + G4UIcommand::ConvertToString(psf_TMin - 1000. * CLHEP::ns) + " ns");
	G4UImanager::GetUIpointer()->ApplyCommand("/g4sipm/persist/tMax " + G4UIcommand::ConvertToString(psf_TMax + 1000. * CLHEP::ns) + " ns");
	// 从文件中读取角度分布
	rnd::RandomDouble *photon_phi = (rnd::RandomDouble *)new rnd::General("resources/psf_phi_distribution.properties");
	rnd::RandomDouble *photon_theta = (rnd::RandomDouble *)new rnd::General("resources/psf_theta_distribution.properties");

	G4UImanager::GetUIpointer()->ApplyCommand(
		boost::str(boost::format("/ps/angle/phiMin %g deg") % (((rnd::General *)photon_phi)->getXMin() / CLHEP::deg)));
	G4UImanager::GetUIpointer()->ApplyCommand(
		boost::str(boost::format("/ps/angle/phiMax %g deg") % (((rnd::General *)photon_phi)->getXMax() / CLHEP::deg)));
	G4UImanager::GetUIpointer()->ApplyCommand(
		boost::str(boost::format("/ps/angle/thetaMin %g deg") % (((rnd::General *)photon_theta)->getXMin() / CLHEP::deg)));
	G4UImanager::GetUIpointer()->ApplyCommand(
		boost::str(boost::format("/ps/angle/thetaMax %g deg") % (((rnd::General *)photon_theta)->getXMax() / CLHEP::deg)));

	// 获取输入的角度分布，如果没有输入，则使用默认的Direction分布
	rnd::ThreeVector *momentum = new rnd::Direction(messenger->getPhiMin(), messenger->getPhiMax(),
													messenger->getThetaMin(), messenger->getThetaMax());

	// Create particles.
	// 生成nParticles个光子
	for (size_t i = 0; i < nParticles; i++)
	{
		particleGun->SetNumberOfParticles(1);
		particleGun->SetParticleDefinition(G4OpticalPhoton::Definition());
		// Dice time.
		// 生成光子的时间
		particleGun->SetParticleTime(time->shoot());
		// Dice energy.
		// 生成光子的能量
		particleGun->SetParticleEnergy(energy->shoot());
		// Dice position.
		// 生成光子的位置
		particleGun->SetParticlePosition(transformToSourcePlane(position->shoot()) + messenger->getPos());
		particleGun->SetParticleMomentumDirection(transformToSourcePlane(momentum->shoot()));
		// Polarization must be computed afterwards.
		// 生成光子的偏振
		particleGun->SetParticlePolarization(
			getPolarizationVector(polar->shoot(), particleGun->GetParticleMomentumDirection()));
		// Fire!
		// 发射光子
		particleGun->GeneratePrimaryVertex(event);
	}
	// Cleanup
	// 释放内存
	delete position;
	delete time;
	delete energy;
	delete polar;
	delete momentum;
}
